GLDS-38 Analysis Using iDEP 0.91

Analysis of GLDS-38 from NASA GeneLab

This R markdown file was auto-generated by the iDEP website Using iDEP 0.91, originally by Steven Xijin.Ge@sdstate.edu

Ge SX, Son EW, Yao R: iDEP: an integrated web application for differential expression and pathway analysis of RNA-Seq data. BMC Bioinformatics 2018, 19(1):534. PMID:30567491

1. Read data

First we set up the working directory to where the files are saved.

 setwd('~/Documents/iDEP/iDEP-GLDS/GLDS38')
 input_biclustMethod <- 'BCCC()' 

R packages used

library(RSQLite,verbose=FALSE)  # for database connection
library(gplots,verbose=FALSE)       # for hierarchical clustering
library(ggplot2,verbose=FALSE)  # graphics
library(e1071,verbose=FALSE)        # computing kurtosis
library(DT,verbose=FALSE)       # for renderDataTable
library(plotly,verbose=FALSE)   # for interactive heatmap
library(reshape2,verbose=FALSE)     # for melt correlation matrix in heatmap

# From Data Read Function
library(edgeR,verbose=FALSE) # count data D.E. 
library(DESeq2,verbose=FALSE) # count data analysis, DEG.DESeq2

# TSNE Plot, tSNEgenePlot 
library(Rtsne,verbose=FALSE)

# PGSA Pathway PGSEA Pathway, PGSEAplot
library(PGSEA,verbose=FALSE)

## Warning: Package 'KEGG.db' is deprecated and will be removed from Bioconductor
##   version 3.12

# DEG.limma
library(limma,verbose=FALSE) # Differential expression
library(statmod,verbose=FALSE)

# enrichment plot
library(dendextend) # customizing tree

# enrich.net2, moduleNetwork
library(igraph)

# Stringdb_geneList, StringDB_GO_enrichmentData, stringDB_network1
# StringDB_network_link
library(STRINGdb,verbose=FALSE)

# gagePathwayData
library(gage,verbose=FALSE) # pathway analysis

#fgseaPathwayData
library(fgsea,verbose=FALSE) # fast GSEA

#ReactomePAPathwayData
library(ReactomePA,verbose=FALSE) # pathway analysis

# KeggImage
library(pathview)

# genomePlot, genomePlotDataPre
library(PREDA,verbose=FALSE)  # showing expression on genome
library(PREDAsampledata,verbose=FALSE) 
library(hgu133plus2.db,verbose=FALSE)

# biclustering
library(biclust,verbose=FALSE)

if(input_biclustMethod == "BCQU()" )
    library(QUBIC,verbose=FALSE) # have trouble installing on Linux
if(input_biclustMethod == "BCUnibic()" )
    library(runibic,verbose=FALSE) # Test biclustMethod dependant qubic runibic

# wgcna
library(WGCNA)
library(flashClust,verbose=FALSE)

 if(file.exists('iDEP_core_functions.R'))
    source('iDEP_core_functions.R') else 
    source('https://raw.githubusercontent.com/iDEP-SDSU/idep/master/shinyapps/idep/iDEP_core_functions.R') 

#Each row of this matrix represents a color scheme;
mycolors = sort(rainbow(20))[c(1,20,10,11,2,19,3,12,4,13,5,14,6,15,7,16,8,17,9,18)]

hmcols <- colorRampPalette(rev(c("#D73027", "#FC8D59", "#FEE090", "#FFFFBF",
"#E0F3F8", "#91BFDB", "#4575B4")))(75)
heatColors = rbind(      greenred(75),     bluered(75),     colorpanel(75,"green","black","magenta"),colorpanel(75,"blue","yellow","red"),hmcols )
rownames(heatColors) = c("Green-Black-Red","Blue-White-Red","Green-Black-Magenta","Blue-Yellow-Red","Blue-white-brown")

We are using the downloaded gene expression file where gene IDs has been converted to Ensembl gene IDs. This is because the ID conversion database is too large to download. You can use your original file if your file uses Ensembl ID, or you do not want to use the pathway files available in iDEP (or it is not available).

 inputFile_ <- 'GLDS38_Expression.csv'
 sampleInfoFile_ <- 'GLDS38_Sampleinfo.csv' 
 gldsMetadataFile_ <- 'GLDS38_Metadata.csv'
 geneInfoFile_ <- 'Arabidopsis_thaliana__athaliana_eg_gene_GeneInfo.csv' #Gene symbols, location etc. 
 geneSetFile_ <- 'Arabidopsis_thaliana__athaliana_eg_gene.db'  # pathway database in SQL; can be GMT format 
 STRING10_speciesFile_ <- 'https://raw.githubusercontent.com/iDEP-SDSU/idep/master/shinyapps/idep/STRING10_species.csv' 

Parameters for reading data

 input_missingValue_ <- 'geneMedian'    #Missing values imputation method
 input_dataFileFormat_ <- 1 #1- read counts, 2 FKPM/RPKM or DNA microarray (Brake up function?)
 input_minCounts_ <- 0.5    #Min counts
 input_NminSamples_ <- 1    #Minimum number of samples 
 input_countsLogStart_ <- 4 #Pseudo count for log CPM
 input_CountsTransform_ <- 1    #Methods for data transformation of counts. 1-EdgeR's logCPM 2-VST, 3-rlog 

readMetadata.out_ <- readMetadata(gldsMetadataFile_)
library(knitr)   #  install if needed. for showing tables with kable
library(kableExtra)
kable( readMetadata.out_ ) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")

	FLT_Rep1	FLT_Rep2	FLT_Rep3	GC_Rep1	GC_Rep2	GC_Rep3	LN2_Rep1	LN2_Rep2	LN2_Rep3	RNAlat_Rep1	RNAlat_Rep2	RNAlat_Rep3
Sample.LongId	Atha.WT.Col.0.sl.FLT.Rep1.G1S1.RNAseq.RNAseq	Atha.WT.Col.0.sl.FLT.Rep2.G1S2.RNAseq.RNAseq	Atha.WT.Col.0.sl.FLT.Rep3.G1S3.RNAseq.RNAseq	Atha.WT.Col.0.sl.GC.Rep1.G2S1.RNAseq.RNAseq	Atha.WT.Col.0.sl.GC.Rep2.G2S2.RNAseq.RNAseq	Atha.WT.Col.0.sl.GC.Rep3.G2S3.RNAseq.RNAseq	Atha.WT.Col.0.sl.LN2.Rep1.n2.1.RNAseq.RNAseq	Atha.WT.Col.0.sl.LN2.Rep2.n2.2.RNAseq.RNAseq	Atha.WT.Col.0.sl.LN2.Rep3.n2.3.RNAseq.RNAseq	Atha.WT.Col.0.sl.RNAlat.Rep1.rl1.RNAseq.RNAseq	Atha.WT.Col.0.sl.RNAlat.Rep2.rl2.RNAseq.RNAseq	Atha.WT.Col.0.sl.RNAlat.Rep3.rl3.RNAseq.RNAseq
Sample.Id	Atha.WT.Col.0.sl.FLT.Rep1.G1S1	Atha.WT.Col.0.sl.FLT.Rep2.G1S2	Atha.WT.Col.0.sl.FLT.Rep3.G1S3	Atha.WT.Col.0.sl.GC.Rep1.G2S1	Atha.WT.Col.0.sl.GC.Rep2.G2S2	Atha.WT.Col.0.sl.GC.Rep3.G2S3	Atha.WT.Col.0.sl.LN2.Rep1.n2.1	Atha.WT.Col.0.sl.LN2.Rep2.n2.2	Atha.WT.Col.0.sl.LN2.Rep3.n2.3	Atha.WT.Col.0.sl.RNAlat.Rep1.rl1	Atha.WT.Col.0.sl.RNAlat.Rep2.rl2	Atha.WT.Col.0.sl.RNAlat.Rep3.rl3
Sample.Name	Atha_WT-Col-0_sl_FLT_Rep1_G1S1	Atha_WT-Col-0_sl_FLT_Rep2_G1S2	Atha_WT-Col-0_sl_FLT_Rep3_G1S3	Atha_WT-Col-0_sl_GC_Rep1_G2S1	Atha_WT-Col-0_sl_GC_Rep2_G2S2	Atha_WT-Col-0_sl_GC_Rep3_G2S3	Atha_WT-Col-0_sl_LN2_Rep1_n2-1	Atha_WT-Col-0_sl_LN2_Rep2_n2-2	Atha_WT-Col-0_sl_LN2_Rep3_n2-3	Atha_WT-Col-0_sl_RNAlat_Rep1_rl1	Atha_WT-Col-0_sl_RNAlat_Rep2_rl2	Atha_WT-Col-0_sl_RNAlat_Rep3_rl3
GLDS	38	38	38	38	38	38	38	38	38	38	38	38
Accession	GLDS-38	GLDS-38	GLDS-38	GLDS-38	GLDS-38	GLDS-38	GLDS-38	GLDS-38	GLDS-38	GLDS-38	GLDS-38	GLDS-38
Hardware	BRIC	BRIC	BRIC	BRIC	BRIC	BRIC	BRIC	BRIC	BRIC	BRIC	BRIC	BRIC
Tissue	Etiolated seedling	Etiolated seedling	Etiolated seedling	Etiolated seedling	Etiolated seedling	Etiolated seedling	Etiolated seedling	Etiolated seedling	Etiolated seedling	Etiolated seedling	Etiolated seedling	Etiolated seedling
Age	8 days	8 days	8 days	8 days	8 days	8 days	8 days	8 days	8 days	8 days	8 days	8 days
Organism	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana	Arabidopsis thaliana
Ecotype	Col-0	Col-0	Col-0	Col-0	Col-0	Col-0	Col-0	Col-0	Col-0	Col-0	Col-0	Col-0
Genotype	WT	WT	WT	WT	WT	WT	WT	WT	WT	WT	WT	WT
Variety	Col-0 WT	Col-0 WT	Col-0 WT	Col-0 WT	Col-0 WT	Col-0 WT	Col-0 WT	Col-0 WT	Col-0 WT	Col-0 WT	Col-0 WT	Col-0 WT
Radiation	Cosmic radiation	Cosmic radiation	Cosmic radiation	Background Earth	Background Earth	Background Earth	Background Earth	Background Earth	Background Earth	Background Earth	Background Earth	Background Earth
Gravity	Microgravity	Microgravity	Microgravity	Terrestrial	Terrestrial	Terrestrial	Terrestrial	Terrestrial	Terrestrial	Terrestrial	Terrestrial	Terrestrial
Developmental	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings	Etiolated 8 day old seedlings
Time.series.or.Concentration.gradient	Single time point	Single time point	Single time point	Single time point	Single time point	Single time point	Single time point	Single time point	Single time point	Single time point	Single time point	Single time point
Light	Dark	Dark	Dark	Dark	Dark	Dark	Dark	Dark	Dark	Dark	Dark	Dark
Assay..RNAseq.	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling	RNAseq Transcription and Proteomic Profiling
Temperature	22-24	22-24	22-24	22-24	22-24	22-24	22-24	22-24	22-24	22-24	22-24	22-24
Treatment.type	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity	Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity
Treatment.intensity	X	X	X	X	X	X	X	X	X	X	X	X
Treament.timing	X	X	X	X	X	X	X	X	X	X	X	X
Preservation.Method.	RNAlater	RNAlater	RNAlater	RNAlater	RNAlater	RNAlater	Liquid Nitrogen	Liquid Nitrogen	Liquid Nitrogen	RNAlater	RNAlater	RNAlater

 readData.out <- readData(inputFile_, input_missingValue=input_missingValue_, input_dataFileFormat = input_dataFileFormat_, input_minCounts=input_minCounts_, input_NminSamples = input_NminSamples_, input_countsLogStart=input_countsLogStart_, input_CountsTransform = input_CountsTransform_) 

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

   kable( head(readData.out$data) ) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%") 

	FLT_Rep1	FLT_Rep2	FLT_Rep3	GC_Rep1	GC_Rep2	GC_Rep3	LN2_Rep1	LN2_Rep2	LN2_Rep3	RNAlat_Rep1	RNAlat_Rep2	RNAlat_Rep3
ATCG00490	19.02600	20.20307	20.47902	19.35267	20.26267	19.81133	17.85944	18.14881	18.28268	19.90399	19.24373	18.83189
AT2G41310	18.59846	19.97237	19.85194	20.18207	18.98753	19.11273	18.29485	18.20631	18.20418	18.69872	18.47971	18.18740
ATCG00020	17.76784	18.79161	19.31738	17.92858	19.00300	18.51242	17.11173	17.35404	17.37928	18.97660	18.68141	18.13132
AT3G21720	17.38043	17.29131	17.40174	18.45645	16.52706	17.41765	16.35671	15.58054	15.69393	14.72961	14.46187	14.90556
AT2G07671	16.80806	17.55228	17.74476	16.89524	17.38578	17.32608	15.95185	15.89930	15.66004	16.40474	16.05359	15.88383
ATCG00280	16.09754	17.02548	17.38910	17.07244	17.48614	17.20404	15.47935	15.34283	15.45231	17.24330	16.73648	16.46234

 readSampleInfo.out <- readSampleInfo(sampleInfoFile_, readData.out) 
 kable( readSampleInfo.out ) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%") 

	Gravity	Preservation.Method.
FLT_Rep1	Microgravity	RNAlater
FLT_Rep2	Microgravity	RNAlater
FLT_Rep3	Microgravity	RNAlater
GC_Rep1	Terrestrial	RNAlater
GC_Rep2	Terrestrial	RNAlater
GC_Rep3	Terrestrial	RNAlater
LN2_Rep1	Terrestrial	Liquid Nitrogen
LN2_Rep2	Terrestrial	Liquid Nitrogen
LN2_Rep3	Terrestrial	Liquid Nitrogen
RNAlat_Rep1	Terrestrial	RNAlater
RNAlat_Rep2	Terrestrial	RNAlater
RNAlat_Rep3	Terrestrial	RNAlater

 input_selectOrg_ ="NEW" 
 input_selectGO_ <- 'GOBP'  #Gene set category 
 input_noIDConversion_ = TRUE  
 allGeneInfo.out_ <- geneInfo(geneInfoFile_) 
 converted.out_ = NULL 
 convertedData.out <- convertedData(converted.out=converted.out_, readData.out=readData.out, input_noIDConversion=input_noIDConversion_)     
 nGenesFilter(readData.out=readData.out, converted.out=converted.out_, convertedData.out=convertedData.out, input_noIDConversion=input_noIDConversion_)  

## [1] "16156 genes in 12 samples. 16117  genes passed filter.\n Original gene IDs used."

 convertedCounts.out_ <- convertedCounts(readData.out=readData.out, converted.out=converted.out_)  # converted counts, just for compatibility 

2. Pre-process

# Read counts per library 
 parDefault = par() 
 par(mar=c(12,4,2,2)) 
 # barplot of total read counts
 x <- readData.out$rawCounts
 groups = as.factor( detectGroups(colnames(x ) ) )
 if(nlevels(groups)<=1 | nlevels(groups) >20 )  
  col1 = 'green'  else
  col1 = rainbow(nlevels(groups))[ groups ]             
         
 barplot( colSums(x)/1e6, 
        col=col1,las=3, main="Total read counts (millions)")  

 readCountsBias(readData.out=readData.out, readSampleInfo.out=readSampleInfo.out)  # detecting bias in sequencing depth 

## [1] 0.02579404
## [1] 0.1547212
## [1] 0.05122098

## [1] "Warning! Sequencing depth bias detected. Total read counts are significantly different among sample groups (p= 2.58e-02 ) based on ANOVA."

 # Box plot 
 x = readData.out$data 
 boxplot(x, las = 2, col=col1,
    ylab='Transformed expression levels',
    main='Distribution of transformed data') 

 #Density plot 
 par(parDefault) 

## Warning in par(parDefault): graphical parameter "cin" cannot be set

## Warning in par(parDefault): graphical parameter "cra" cannot be set

## Warning in par(parDefault): graphical parameter "csi" cannot be set

## Warning in par(parDefault): graphical parameter "cxy" cannot be set

## Warning in par(parDefault): graphical parameter "din" cannot be set

## Warning in par(parDefault): graphical parameter "page" cannot be set

 densityPlot(readData.out = readData.out, mycolors)       

 # Scatter plot of the first two samples 
 plot(x[,1:2],xlab=colnames(x)[1],ylab=colnames(x)[2], 
    main='Scatter plot of first two samples') 

 ####plot gene or gene family
 input_selectOrg_ ="BestMatch" 
 input_geneSearch_ <- 'HOXA'    #Gene ID for searching 
 genePlot(allGeneInfo.out_, convertedData.out=convertedData.out, input_selectOrg=input_selectOrg_, input_geneSearch=input_geneSearch_)  

## NULL

 input_useSD_ <- 'FALSE'    #Use standard deviation instead of standard error in error bar? 
 geneBarPlotError(allGeneInfo.out_, input_selectOrg_, input_geneSearch_, input_useSD_)       

## NULL

3. Heatmap

 # hierarchical clustering tree
 x <- readData.out$data
 maxGene <- apply(x,1,max)
 # remove bottom 25% lowly expressed genes, which inflate the PPC
 x <- x[which(maxGene > quantile(maxGene)[1] ) ,] 
 plot(as.dendrogram(hclust2( dist2(t(x)))), ylab="1 - Pearson C.C.", type = "rectangle") 

 #Correlation matrix
 input_labelPCC_ <- TRUE    #Show correlation coefficient? 
 correlationMatrix(input_labelPCC = input_labelPCC_) 

 # Parameters for heatmap
 input_nGenes_ <- 1000  #Top genes for heatmap
 input_geneCentering_ <- TRUE   #centering genes ?
 input_sampleCentering_ <- FALSE    #Center by sample?
 input_geneNormalize_ <- FALSE  #Normalize by gene?
 input_sampleNormalize_ <- FALSE    #Normalize by sample?
 input_noSampleClustering_ <- FALSE #Use original sample order
 input_heatmapCutoff_ <- 4  #Remove outliers beyond number of SDs 
 input_distFunctions_ <- 1  #which distant funciton to use
 input_hclustFunctions_ <- 1    #Linkage type
 input_heatColors1_ <- 1    #Colors
 input_selectFactorsHeatmap_ <- 'Gravity'   #Sample coloring factors 
 png('heatmap.png', width = 10, height = 15, units = 'in', res = 300) 
 staticHeatmap(heatColors, input_nGenes_, input_geneCentering_, input_sampleCentering_, input_geneNormalize_, input_sampleNormalize_, input_noSampleClustering=input_noSampleClustering_, input_heatmapCutoff=input_heatmapCutoff_, input_distFunctions=input_distFunctions_, input_hclustFunctions=input_hclustFunctions_, input_heatColors1=input_heatColors1_, input_selectFactorsHeatmap=input_selectFactorsHeatmap_) 

 dev.off()  

## png 
##   2

[heatmap] (heatmap.png)

 heatmapPlotly(heatColors, allGeneInfo.out=allGeneInfo.out_, input_geneCentering=input_geneCentering_, input_sampleCentering = input_sampleCentering_, input_geneNormalize = input_geneNormalize_, input_sampleNormalize = input_sampleNormalize_, input_heatColors1=input_heatColors1_) # interactive heatmap using Plotly 

4. K-means clustering

 input_nGenesKNN_ <- 2000   #Number of genes fro k-Means
 input_nClusters_ <- 4  #Number of clusters 
 maxGeneClustering_ = 12000
 input_kmeansNormalization_ <- 'geneMean'   #Normalization
 input_KmeansReRun_ <- 0    #Random seed 

 distributionSD(input_nGenesKNN_)  #Distribution of standard deviations 

 KmeansNclusters(input_nGenesKNN_)  #Number of clusters 

 Kmeans.out = Kmeans(maxGeneClustering=maxGeneClustering_, input_nGenesKNN=input_nGenesKNN_, input_nClusters=input_nClusters_, input_kmeansNormalization=input_kmeansNormalization_, input_KmeansReRun=input_KmeansReRun_)   #Running K-means 
 KmeansHeatmap(.mycolors = mycolors, .heatColors = heatColors, .input_heatColors1=input_heatColors1_)   #Heatmap for k-Means 

 #Read gene sets for enrichment analysis 
 sqlite_  <- dbDriver('SQLite') #sqlite changed to default in most functions
 input_selectGO3_ <- 'GOBP' #Gene set category
 input_minSetSize_ <- 15    #Min gene set size
 input_maxSetSize_ <- 2000  #Max gene set size 
 GeneSets.out <-readGeneSets( geneSetFile_,
    convertedData.out, input_selectGO3_,input_selectOrg_,
    c(input_minSetSize_, input_maxSetSize_)  )  
 # Alternatively, users can use their own GMT files by
 #GeneSets.out <- readGMTRobust('somefile.GMT')  
 results <- KmeansGO(input_nClusters=input_nClusters_)  #Enrichment analysis for k-Means clusters   
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%") 

Cluster	adj.Pval	Genes	Pathways
A	7.99e-122	219	Organonitrogen compound biosynthetic process
	2.84e-118	155	Translation
	4.26e-118	155	Peptide biosynthetic process
	1.81e-116	159	Amide biosynthetic process
	3.07e-113	156	Peptide metabolic process
	1.22e-106	160	Cellular amide metabolic process
	5.05e-51	157	Response to abiotic stimulus
	4.90e-42	57	Photosynthesis
	5.59e-35	50	Response to cytokinin
	7.30e-34	63	Generation of precursor metabolites and energy
B	9.12e-48	81	Cell wall organization or biogenesis
	9.12e-48	73	Cell wall organization
	2.49e-47	74	External encapsulating structure organization
	1.20e-40	128	Response to abiotic stimulus
	2.58e-36	74	Drug metabolic process
	8.64e-36	111	Small molecule metabolic process
	2.30e-34	82	Carbohydrate metabolic process
	2.27e-32	77	Response to inorganic substance
	6.60e-28	42	Plant-type cell wall organization or biogenesis
	1.13e-27	51	Polysaccharide metabolic process
C	3.98e-50	79	Response to external stimulus
	1.34e-42	63	Response to external biotic stimulus
	1.34e-42	63	Response to other organism
	2.33e-42	63	Response to biotic stimulus
	4.56e-40	42	Secondary metabolic process
	1.25e-38	63	Defense response
	2.03e-38	67	Multi-organism process
	2.35e-35	78	Response to abiotic stimulus
	4.27e-34	26	Indole-containing compound metabolic process
	1.99e-33	49	Defense response to other organism
D	1.20e-50	151	Response to abiotic stimulus
	1.95e-47	127	Response to oxygen-containing compound
	2.73e-46	98	Response to inorganic substance
	1.10e-39	101	Response to acid chemical
	2.11e-39	113	Cellular response to chemical stimulus
	6.77e-37	126	Response to organic substance
	1.17e-31	107	Response to hormone
	4.57e-31	107	Response to endogenous stimulus
	3.63e-27	75	Response to external biotic stimulus
	3.63e-27	75	Response to other organism

 input_seedTSNE_ <- 0   #Random seed for t-SNE
 input_colorGenes_ <- TRUE  #Color genes in t-SNE plot? 
 tSNEgenePlot(input_seedTSNE=input_seedTSNE_, input_colorGenes=input_colorGenes_)  #Plot genes using t-SNE 

5. PCA and beyond

 input_selectFactors_ <- 'Gravity'  #Factor coded by color
 input_selectFactors2_ <- 'Preservation.Method.'    #Factor coded by shape
 input_tsneSeed2_ <- 0  #Random seed for t-SNE 
 #PCA, MDS and t-SNE plots
 PCAplot(input_selectFactors_, input_selectFactors2_)  

 MDSplot(input_selectFactors_, input_selectFactors2_) 

 tSNEplot(input_selectFactors_, input_selectFactors2_, input_tsneSeed2_)  

 #Read gene sets for pathway analysis using PGSEA on principal components 
 input_selectGO6_ <- 'GOBP' 
 GeneSets.out_ <-readGeneSets( geneSetFile_,
    convertedData.out, input_selectGO6_,input_selectOrg_,
    c(input_minSetSize_, input_maxSetSize_)  )  
 PCApathway(convertedData.out=convertedData.out, GeneSets.out=GeneSets.out_) # Run PGSEA analysis 

 cat( PCA2factor(readData.out, readSampleInfo.out=readSampleInfo.out) )   #The correlation between PCs with factors 

## 
##  Correlation between Principal Components (PCs) with factors
## PC1 is correlated with Gravity (p=3.91e-02).

6. DEG1

 input_CountsDEGMethod_ <- 3    #DESeq2= 3,limma-voom=2,limma-trend=1 
 input_limmaPval_ <- 0.1    #FDR cutoff
 input_limmaFC_ <- 2    #Fold-change cutoff
 input_selectModelComprions_ <- 'Gravity: Microgravity vs. Terrestrial' #Selected comparisons
 input_selectFactorsModel_ <- 'Gravity' #Selected comparisons
 input_selectInteractions_ <- NULL  #Selected comparisons
 input_selectBlockFactorsModel_ <- NULL     #Selected comparisons
 factorReferenceLevels.out_ <- c('Gravity:Terrestrial') 

 limma.out <- limma(input_dataFileFormat = input_dataFileFormat_, input_countsLogStart = input_countsLogStart_, convertedCounts.out=convertedCounts.out_, input_CountsDEGMethod=input_CountsDEGMethod_, input_limmaPval=input_limmaPval_, input_limmaFC=input_limmaFC_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel=input_selectFactorsModel_, input_selectInteractions=input_selectInteractions_, input_selectBlockFactorsModel=input_selectBlockFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

 DEG.data.out <- DEG.data(allGeneInfo.out_)
 limma.out$comparisons 

## [1] "Microgravity-Terrestrial"

 input_selectComparisonsVenn_ = limma.out$comparisons[1:3] # use first three comparisons
 input_UpDownRegulated_ <- FALSE    #Split up and down regulated genes 
 vennPlot(input_selectComparisonsVenn_, input_UpDownRegulated_) # Venn diagram 

  sigGeneStats() # number of DEGs as figure 

  sigGeneStatsTable() # number of DEGs as table 

##                                       Comparisons  Up Down
## Microgravity-Terrestrial Microgravity-Terrestrial 353 1071

7. DEG2

 input_selectContrast_ <- 'Microgravity-Terrestrial'    #Selected comparisons 
 selectedHeatmap.data.out <- selectedHeatmap.data(.converted.out=converted.out_, .readData.out=readData.out, .input_noIDConversion=input_noIDConversion_, input_dataFileFormat=input_dataFileFormat_, input_CountsDEGMethod=input_CountsDEGMethod_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel=input_selectFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_, input_selectContrast=input_selectContrast_)
 selectedHeatmap(.mycolors=mycolors,.heatColors=heatColors)   # heatmap for DEGs in selected comparison

 # Save gene lists and data into files
 write.csv( selectedHeatmap.data(.converted.out=converted.out_, .readData.out=readData.out, .input_noIDConversion=input_noIDConversion_, input_dataFileFormat = input_dataFileFormat_, input_CountsDEGMethod=input_CountsDEGMethod_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel=input_selectFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_, input_selectContrast=input_selectContrast_)$genes, 'heatmap.data.csv') 
 write.csv(DEG.data(allGeneInfo.out_),'DEG.data.csv' )
 write(AllGeneListsGMT() ,'AllGeneListsGMT.gmt')

 input_selectGO2_ <- 'GOBP' #Gene set category 
 geneListData.out <- geneListData(allGeneInfo.out_, input_selectGO2_, input_selectOrg_, input_limmaPval_, input_limmaFC_, input_selectContrast_)  
 volcanoPlot(input_limmaPval=input_limmaPval_, input_limmaFC=input_limmaFC_, input_selectContrast=input_selectContrast_)  

  scatterPlot(input_dataFileFormat=input_dataFileFormat_, input_CountsDEGMethod=input_CountsDEGMethod_, input_limmaPval=input_limmaPval_, input_limmaFC=input_limmaFC_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel=input_selectFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_, input_selectContrast_)  

  MAplot(.converted.out=converted.out_, .readData.out=readData.out, .input_noIDConversion=input_noIDConversion_, input_dataFileFormat=input_dataFileFormat_, input_CountsDEGMethod=input_CountsDEGMethod_, input_limmaPval=input_limmaPval_, input_limmaFC=input_limmaFC_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel_, factorReferenceLevels.out_, input_selectContrast_)  

  geneListGOTable.out <- geneListGOTable()  
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile_,
    convertedData.out, input_selectGO2_,input_selectOrg_,
    c(input_minSetSize_, input_maxSetSize_)  ) 
 input_removeRedudantSets_ <- TRUE  #Remove highly redundant gene sets? 
 results <- geneListGO(input_removeRedudantSets_)  #Enrichment analysis
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%") 

Direction	adj.Pval	nGenes	Pathways
Down regulated	4.1e-37	143	Response to external stimulus
	4.6e-36	117	Response to biotic stimulus
	7.5e-35	114	Response to external biotic stimulus
	7.5e-35	114	Response to other organism
	1.1e-34	124	Defense response
	2.3e-33	134	Multi-organism process
	1.7e-28	60	Secondary metabolic process
	1.1e-27	79	Response to drug
	1.5e-26	137	Response to oxygen-containing compound
	2.8e-26	156	Response to organic substance
Up regulated	1.9e-18	29	Photosynthesis
	5.7e-15	61	Organonitrogen compound biosynthetic process
	7.5e-15	41	Peptide metabolic process
	5.2e-14	38	Translation
	5.2e-14	38	Peptide biosynthetic process
	1.6e-13	42	Cellular amide metabolic process
	1.6e-13	39	Amide biosynthetic process
	3.1e-13	19	Pigment biosynthetic process
	1.1e-12	20	Pigment metabolic process
	1.3e-11	23	Plastid organization

STRING-db API access. We need to find the taxonomy id of your species, this used by STRING. First we try to guess the ID based on iDEP’s database. Users can also skip this step and assign NCBI taxonomy id directly by findTaxonomyID.out = 10090 # mouse 10090, human 9606 etc.

 STRING10_species_ = read.csv(STRING10_speciesFile_)  
 ix = grep('Arabidopsis thaliana', STRING10_species_$official_name ) 
 findTaxonomyID.out <- STRING10_species_[ix,1] # find taxonomyID
 findTaxonomyID.out  

## [1] 3702

Enrichment analysis using STRING

  STRINGdb_geneList.out <- STRINGdb_geneList() #convert gene lists

## Warning:  we couldn't map to STRING 0% of your identifiers

 input_STRINGdbGO_ <- 'Process' #'Process', 'Component', 'Function', 'KEGG', 'Pfam', 'InterPro' 
 results <- stringDB_GO_enrichmentData(input_STRINGdbGO_)  # enrichment using STRING    
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%") 

“No significant enrichment found.”	adj.Pval
No significant enrichment found.	NULL

PPI network retrieval and analysis

 input_nGenesPPI_ <- 100    #Number of top genes for PPI retrieval and analysis 
 stringDB_network1(1, input_nGenesPPI_) #Show PPI network 

Generating interactive PPI

 write(stringDB_network_link(input_nGenesPPI_), 'PPI_results.html') # write results to html file 

## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")

## Warning:  we couldn't map to STRING 0% of your identifiers

## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")

## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")

 browseURL('PPI_results.html') # open in browser 

8. Pathway analysis

 input_selectContrast1_ <- 'Microgravity-Terrestrial'   #select Comparison 
 #input_selectContrast1 = limma.out$comparisons[3] # manually set
 input_selectGO_ <- 'GOBP'  #Gene set category 
 #input_selectGO='custom' # if custom gmt file
 input_minSetSize_ <- 15    #Min size for gene set
 input_maxSetSize_ <- 2000  #Max size for gene set 
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile_,
    convertedData.out, input_selectGO_,input_selectOrg_,
    c(input_minSetSize_, input_maxSetSize_)  ) 
 input_pathwayPvalCutoff_ <- 0.2    #FDR cutoff
 input_nPathwayShow_ <- 30  #Top pathways to show
 input_absoluteFold_ <- FALSE   #Use absolute values of fold-change?
 input_GenePvalCutoff_ <- 1 #FDR to remove genes 

 input_pathwayMethod_ = 1  # 1  GAGE
 gagePathwayData.out <- gagePathwayData(input_minSetSize=input_minSetSize_, input_maxSetSize=input_maxSetSize_, input_selectContrast1=input_selectContrast1_, input_pathwayPvalCutoff_, input_nPathwayShow_, input_absoluteFold_, input_GenePvalCutoff_)  # pathway analysis using GAGE  
   
 results <- gagePathwayData.out  #Enrichment analysis for k-Means clusters  
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%") 

Direction	GAGE analysis: Microgravity vs Terrestrial	statistic	Genes	adj.Pval
Down	Response to chitin	-7.7014	104	1.9e-10
	Secondary metabolic process	-7.4839	265	1.9e-10
	Response to drug	-7.3469	481	1.9e-10
	Response to organonitrogen compound	-6.3515	212	1.3e-07
	Secondary metabolite biosynthetic process	-6.3149	118	2.4e-07
	Response to nitrogen compound	-6.2461	271	1.6e-07
Up	Photosynthesis	10.7757	223	3.7e-21
	Ribosome biogenesis	10.2165	342	6.5e-20
	Ribonucleoprotein complex biogenesis	10.0934	437	6.5e-20
	NcRNA metabolic process	9.4896	425	8.0e-18
	Plastid organization	9.039	257	6.2e-16
	NcRNA processing	8.7867	357	2.8e-15
	RRNA processing	8.2838	239	2.6e-13
	RRNA metabolic process	8.2747	244	2.6e-13
	Photosynthesis, light reaction	7.9994	119	7.6e-12
	Chloroplast organization	7.4927	198	4.2e-11
	Generation of precursor metabolites and energy	7.2638	391	7.9e-11
	RNA modification	7.1316	321	2.8e-10
	Thylakoid membrane organization	6.8664	46	1.5e-07
	Tetrapyrrole biosynthetic process	6.0669	70	7.8e-07
	Ribosome assembly	6.0176	76	1.3e-06
	Ribosomal large subunit biogenesis	5.9801	99	1.0e-06
	Tetrapyrrole metabolic process	5.964	93	7.8e-07
	Porphyrin-containing compound biosynthetic process	5.9033	67	1.3e-06
	Porphyrin-containing compound metabolic process	5.8821	92	1.1e-06
	Protein transmembrane transport	5.7744	112	1.4e-06
	Ribonucleoprotein complex subunit organization	5.7385	181	1.3e-06
	Nucleoside monophosphate metabolic process	5.6492	193	1.6e-06
	Ribonucleoprotein complex assembly	5.6323	173	1.8e-06
	TRNA metabolic process	5.6148	170	1.8e-06

 pathwayListData.out_ = pathwayListData(allGeneInfo.out_, input_selectOrg_, input_selectGO_, input_pathwayMethod_) 
 enrichmentPlot(pathwayListData.out_, 25  ) 

  enrichmentNetwork(pathwayListData.out_ )  

  enrichmentNetworkPlotly(pathwayListData.out_) 

## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

 input_pathwayMethod_ = 3  # 1  fgsea 
 fgseaPathwayData.out <- fgseaPathwayData(input_minSetSize=input_minSetSize_, input_maxSetSize=input_maxSetSize_, input_selectContrast1=input_selectContrast1_, input_pathwayPvalCutoff=input_pathwayPvalCutoff_, input_absoluteFold=input_absoluteFold_, input_nPathwayShow=input_nPathwayShow_, input_GenePvalCutoff=input_GenePvalCutoff_) #Pathway analysis using fgsea 

## Warning in fgsea(pathways = gmt, stats = fold, minSize = input_minSetSize, :
## You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To
## run fgseaMultilevel, you need to remove the nperm argument in the fgsea function
## call.

## Warning in fgseaSimple(...): There were 3 pathways for which P-values were not
## calculated properly due to unbalanced gene-level statistic values

 results <- fgseaPathwayData.out  #Enrichment analysis for k-Means clusters 
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%") 

Direction	GSEA analysis: Microgravity vs Terrestrial	NES	Genes	adj.Pval
Up	Photosynthesis	3.2975	223	1.2e-02
	Photosynthesis, light reaction	3.0161	119	6.3e-03
	Thylakoid membrane organization	2.9239	46	4.1e-03
	Translation	2.9137	619	1.1e-01
	Peptide biosynthetic process	2.9053	622	1.1e-01
	Plastid organization	2.898	257	1.4e-02
	Amide biosynthetic process	2.8546	683	1.8e-01
	Ribosome biogenesis	2.8254	342	2.3e-02
	Photosynthetic electron transport chain	2.8113	46	4.1e-03
	Chloroplast organization	2.7408	198	9.8e-03
	Ribosome assembly	2.7189	76	4.8e-03
	Tetrapyrrole biosynthetic process	2.6943	70	4.7e-03
	RRNA processing	2.6919	239	1.3e-02
	Plastid membrane organization	2.6845	49	4.2e-03
	RRNA metabolic process	2.6757	244	1.3e-02
	Chlorophyll biosynthetic process	2.6734	58	4.3e-03
	Protein localization to chloroplast	2.6602	45	4.1e-03
	Porphyrin-containing compound biosynthetic process	2.6383	67	4.6e-03
	Ribonucleoprotein complex biogenesis	2.6344	437	3.9e-02
	Ribosomal large subunit biogenesis	2.6273	99	5.6e-03
	Tetrapyrrole metabolic process	2.6198	93	5.4e-03
	Protein targeting to chloroplast	2.6126	43	4.1e-03
	Establishment of protein localization to chloroplast	2.6126	43	4.1e-03
	Chlorophyll metabolic process	2.6057	81	4.9e-03
	Porphyrin-containing compound metabolic process	2.5984	92	5.4e-03
	Peptide metabolic process	2.5698	679	1.7e-01
	Chloroplast RNA processing	2.5648	19	3.9e-03
	Pigment biosynthetic process	2.5423	129	6.8e-03
	NcRNA processing	2.499	357	2.4e-02
	Chloroplast rRNA processing	2.4985	18	3.9e-03

  pathwayListData.out_ = pathwayListData(allGeneInfo.out_, input_selectOrg=input_selectOrg_, input_selectGO=input_selectGO_, input_pathwayMethod_) 
  enrichmentPlot(pathwayListData.out_, 25, mycolors) 

  enrichmentNetwork(pathwayListData.out_ )  

  enrichmentNetworkPlotly(pathwayListData.out_) 

   PGSEAplot(input_selectOrg = input_selectOrg_, input_dataFileFormat=input_dataFileFormat_, input_selectGO=input_selectGO_, input_minSetSize=input_minSetSize_, input_maxSetSize=input_maxSetSize_, input_CountsDEGMethod=input_CountsDEGMethod_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_, input_selectContrast1=input_selectContrast1_, input_pathwayPvalCutoff=input_pathwayPvalCutoff_, input_nPathwayShow_) # pathway analysis using PGSEA 

## 
## Computing P values using ANOVA

9. Chromosome

 input_selectContrast2_ <- 'Microgravity-Terrestrial'   #select Comparison 
 #input_selectContrast2 = limma.out$comparisons[3] # manually set
 input_limmaPvalViz_ <- 0.1 #FDR to filter genes
 input_limmaFCViz_ <- 2 #FDR to filter genes 
 genomePlotly(allGeneInfo.out_, input_selectContrast2_, input_limmaPvalViz_, input_limmaFCViz_) # shows fold-changes on the genome 

## Warning in eval(quote(list(...)), env): NAs introduced by coercion

## Warning in genomePlotly(allGeneInfo.out_, input_selectContrast2_,
## input_limmaPvalViz_, : NAs introduced by coercion

10. Biclustering

 input_nGenesBiclust_ <- 1000   #Top genes for biclustering
 input_biclustMethod_ <- 'BCCC()'   #Method: 'BCCC', 'QUBIC', 'runibic' ... 
 biclustering.out = biclustering(input_nGenesBiclust_, input_biclustMethod_ )  # run analysis

 input_selectBicluster_ <- 1    #select a cluster 
 biclustHeatmap(heatColors, input_heatColors1=input_heatColors1_, input_selectBicluster=input_selectBicluster_)   # heatmap for selected cluster 

 input_selectGO4_ <- 'GOBP' #Gene set category 
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile_,
    convertedData.out, input_selectGO4_, input_selectOrg_,
    c(input_minSetSize_, input_maxSetSize_)  )  
 results <- geneListBclustGO(2, input_selectBicluster=input_selectBicluster_)  #Enrichment analysis for k-Means clusters    
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%") 

adj.Pval	Genes	Pathways
2.4e-134	308	Response to abiotic stimulus
5.9e-86	170	Response to inorganic substance
1.1e-66	216	Response to organic substance
5.7e-64	182	Oxidation-reduction process
1.3e-63	197	Organonitrogen compound biosynthetic process
2.2e-63	173	Response to external stimulus
4.5e-62	144	Response to external biotic stimulus
4.5e-62	144	Response to other organism
2.6e-61	144	Response to biotic stimulus
3.1e-61	105	Response to metal ion

11. Co-expression network

 input_mySoftPower_ <- 5    #SoftPower to cutoff
 input_nGenesNetwork_ <- 1000   #Number of top genes
 input_minModuleSize_ <- 20 #Module size minimum 
 wgcna.out = wgcna(2000, input_mySoftPower=input_mySoftPower_, input_nGenesNetwork=input_nGenesNetwork_, input_minModuleSize=input_minModuleSize_)   # run WGCNA  

## Warning: executing %dopar% sequentially: no parallel backend registered

##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.6880  3.020          0.962  418.00    426.00  557.0
## 2      2   0.5030  1.310          0.949  240.00    241.00  388.0
## 3      3   0.2260  0.549          0.896  159.00    155.00  294.0
## 4      4   0.0337  0.163          0.868  113.00    108.00  233.0
## 5      5   0.0639 -0.187          0.865   85.00     79.00  191.0
## 6      6   0.2370 -0.350          0.875   66.30     60.10  159.0
## 7      7   0.4190 -0.473          0.908   53.20     46.40  135.0
## 8      8   0.5660 -0.576          0.946   43.70     37.00  117.0
## 9      9   0.6550 -0.661          0.944   36.50     30.10  103.0
## 10    10   0.7000 -0.706          0.948   30.90     24.80   91.1
## 11    12   0.7590 -0.822          0.928   23.00     17.50   73.4
## 12    14   0.7970 -0.918          0.922   17.70     12.80   60.7
## 13    16   0.8360 -0.968          0.932   14.00      9.68   51.3
## 14    18   0.8470 -1.020          0.929   11.40      7.55   44.1
## 15    20   0.8340 -1.050          0.905    9.39      5.96   38.5
## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.

 softPower(wgcna.out)  # soft power curve 

  modulePlot(wgcna.out)  # plot modules  

  listWGCNA.Modules.out = listWGCNA.Modules() #modules

 input_selectGO5_ <- 'GOBP' #Gene set category 
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile_,
    convertedData.out, input_selectGO5_,input_selectOrg_,
    c(input_minSetSize_, input_maxSetSize_)  ) 
 input_selectWGCNA.Module_ <- 'Entire network'  #Select a module
 input_topGenesNetwork_ <- 10   #SoftPower to cutoff
 input_edgeThreshold_ <- 0.4    #Number of top genes 
 moduleNetwork(input_noIDConversion=input_noIDConversion_, allGeneInfo.out=allGeneInfo.out_, input_selectOrg_, input_selectWGCNA.Module_, input_topGenesNetwork_, input_edgeThreshold_, input_selectGO5_)   # show network of top genes in selected module

##  softConnectivity: FYI: connecitivty of genes with less than 4 valid samples will be returned as NA.
##  ..calculating connectivities..

 input_removeRedudantSets_ <- TRUE  #Remove redundant gene sets 
 results <- networkModuleGO(2, input_selectWGCNA.Module_)  #Enrichment analysis of selected module
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%") 

adj.Pval	Genes	Pathways
2.4e-134	308	Response to abiotic stimulus
5.9e-86	170	Response to inorganic substance
1.1e-66	216	Response to organic substance
5.7e-64	182	Oxidation-reduction process
1.3e-63	197	Organonitrogen compound biosynthetic process
2.2e-63	173	Response to external stimulus
4.5e-62	144	Response to external biotic stimulus
4.5e-62	144	Response to other organism
2.6e-61	144	Response to biotic stimulus
3.1e-61	105	Response to metal ion